import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)
if root == '/Users/mo2016':
    modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
    modelling_home = '/Volumes/mo2016/home/Documents/modelling'
    modelling_local = root + '/Documents/modelling'
    import matplotlib as mpl
    mpl.use('tkagg')

if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
        modelling_ephemeral = root + '/ephemeral/Documents/modelling'
        modelling_home = root  + '/home/Documents/modelling'
        modelling_local = modelling_home

modulepath = modelling_local + '/3954/modules/new_CN'
print(root)
sys.path.append(modulepath)
from plotting_numerical import *
from tqdm import tqdm
import pickle
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import matplotlib.animation as animation
import matplotlib as mpl
#############################
#Opening list with parID's

lhs=True
if lhs==True:
    folder = 'fullcircuit_5716gaussian/1M_turingI'
    variant=0
    details = '1M_turingI'
else:
    var=float(sys.argv[1])
    # var=0.23
    folder = 'fullcircuit_5716gaussian/var%s'%var
    variant='5716gaussian'
    details = 'var%r'%var

circuit_n=2
shape='ca'
mechanism = 'fullcircuit'
dimension='2D'


L=10; x_gridpoints =15; J = L*x_gridpoints


T =120; t_gridpoints = 10; N = T*t_gridpoints
data_path = modelling_home + '/3954/numerical_confocal/results/simulation/ca/%s'%(folder)
if lhs==True:    
    if folder == 'fullcircuit_5716gaussian/1M_turingI':
        print('yes')
        general_filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,'1MturingI', shape,mechanism,L,J,T,N)
    else: 
        general_filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
    general_filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)
# metric=str(sys.argv[1])

metric='ps_min'
parID_dict = pickle.load( open(modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/%s_dict_%s.pkl"%(metric,general_filename), "rb" ) )
print(modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/%s_dict_%s.pkl"%(metric,general_filename))
len(parID_dict)
parID_list = []
entropy_list = []
parID_dict = dict(sorted(parID_dict.items(), key=lambda item: item[1])) #important to sort out dictionary by HKS values
for key in parID_dict:
    parID_list.append(key)
    entropy_list.append(parID_dict[key])
# print(parID_list, entropy_list)

parID_list = [int(i) for i in parID_list] #turn string list into integer list
print(parID_list)
# k=20

num=len(parID_list)
n_col = int(np.sqrt(num))
# for i in range(0,k):
# n_col = 10
#     row = np.where(fit==i)[0]
#     print(row) # row in Z for elements of cluster i
#     num = row.shape[0]       #  number of elements for each cluster

n_row = np.floor(num/n_col)+1    # number of rows in the figure of the cluster
#     print("cluster "+str(i))
#     print(str(num)+" elements")

fig = plt.figure(figsize=(n_col/10+2, n_row/10+2))
dx = float(L)/float(J-1)
grid = np.array([j*dx for j in range(J)])
for count,parID in tqdm(enumerate(parID_list[:10]),disable=False):
    ax=plt.subplot(n_row,n_col, count+1)
    #rgb_timeseries=timeseries_unstacked_list[row[n]] # Read the numpy matrix with images in the rows
    # par_ID = parID_list[row[n]]
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
    final_concentration = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
    rgb = plot_redgreenblue_contrast(final_concentration,L,mechanism,shape,filename,modelling_ephemeral,parID=parID,scale_factor=x_gridpoints,save_figure='LargeImage')
    # rgb_timeseries=timeseries_unstacked # Read the numpy matrix with images in the rows
    # ax.set_title(parID,size=0.1)
    ax.set(yticks=[],xticks=[],yticklabels=[],xticklabels=[])
    ax.imshow(rgb.astype('uint8'), origin= 'lower')
    ax.set_ylabel('%s-%s'%(parID,entropy_list[count]),size= 1,c='y', labelpad=0.35)
# plt.title('1M numerical search 0-%r'%num)
# plt.savefig(modelling_home + '/3954/numerical_confocal/results/entropy/LargeImages/%s_%s.png'%(metric,general_filename), dpi=2000)

plt.show()
# plt.savefig('h.png',dpi=2000)

print('gh')
# plt.clf()
# plt.close(fig)

Running cells with 'Python 3.7.3 64-bit' requires ipykernel package.

Run the following command to install 'ipykernel' into the Python environment. 

Command: '/usr/bin/python3 -m pip install ipykernel -U --user --force-reinstall'

import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)
if root == '/Users/mo2016':
    modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
    modelling_home = '/Volumes/mo2016/home/Documents/modelling'
    modelling_local = root + '/Documents/modelling'
    import matplotlib as mpl
    mpl.use('tkagg')

if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
        modelling_ephemeral = root + '/ephemeral/Documents/modelling'
        modelling_home = root  + '/home/Documents/modelling'
        modelling_local = modelling_home

modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)
from plotting_numerical import plot_redgreen_contrast
# from sendmail import *

from tqdm import tqdm
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################
# %matplotlib inline

def fourier2d(final_concentration):
    # print('Normalised PowerSpectrum of original')
    pixels=final_concentration[5]
    pixels[pixels == 0] = 0.0000000001
    pixels/=np.sum(pixels) #line added
    #fourier transform
    f = np.fft.ifftshift(pixels)
    f = np.fft.fft2(f)
    f = np.fft.fftshift(f)
    if plot==True:
        plt.subplot(142)
        plt.imshow(np.real(f))
        plt.colorbar()
        plt.subplot(143)
        plt.imshow(np.real(f))
        plt.colorbar()
        # max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
    PowerSpectrum = np.real(f)**2 +  np.imag(f)**2
    # PowerSpectrum[max_coordinates]=0
    if plot==True:
        plt.subplot(144)
        plt.imshow(PowerSpectrum)
        plt.colorbar()
    if np.sum(PowerSpectrum)!=0:
        PowerSpectrum/=np.sum(PowerSpectrum)
    # print('sum',np.sum(PowerSpectrum))
    # print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))
    # plt.subplot(132)
    # plt.imshow(PowerSpectrum)#,norm=LogNorm())
    # plt.colorbar()
    # binwidth = 0.00001
    PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
    # print(np.sum(PowerSpectrum_vector))
    # plt.subplot(133)
    # hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
    # plt.yscale('log')
    # plt.show()
    return PowerSpectrum_vector

def entropy(histogram):
    x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
    return x

lhs=True
if lhs==True:
    folder = 'fullcircuit_5716gaussian/1M_turingI'
    variant=0
else:
    var=float(sys.argv[1])
    # var=0.23
    folder = 'fullcircuit_5716gaussian/var%s'%var
    variant='5716gaussian'

circuit_n=2
shape='square'
mechanism = 'fullcircuit'


L=5; x_gridpoints =10; J = L*x_gridpoints


T =2000; t_gridpoints = 10; N = T*t_gridpoints
data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )
parID_ps = {}
plot=False
# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
# for parID in tqdm(parID_list[:20], disable=False):

for parID in tqdm([805,686,472,252,688], disable=True):
    # print('parID',parID)
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
    final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
    if not np.isnan(final_concentration0).any():
        final_concentration = np.round(final_concentration0,4)
        if plot==True:
            plt.figure(figsize=[14,2])
            plt.subplot(141)
            plt.imshow(final_concentration[5])
            plt.colorbar()
        hist_PowerSpectrum  = fourier2d(final_concentration)
        ps = entropy(hist_PowerSpectrum)
        parID_ps[parID] = ps
    # print(final_concentration0)

if lhs==True:
    filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
    filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)

pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )
# if root=='/rds/general/user/mo2016':
#     sendmail('general_filename')

Running cells with 'Python 3.7.3 64-bit' requires ipykernel package.

Run the following command to install 'ipykernel' into the Python environment. 

Command: '/usr/bin/python3 -m pip install ipykernel -U --user --force-reinstall'

Python 3.8.5 (default, Sep 4 2020, 07:30:14)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.19.0 -- An enhanced Interactive Python. Type '?' for help.

import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)
if root == '/Users/mo2016':
    modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
    modelling_home = '/Volumes/mo2016/home/Documents/modelling'
    modelling_local = root + '/Documents/modelling'
    import matplotlib as mpl
    mpl.use('tkagg')

if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
        modelling_ephemeral = root + '/ephemeral/Documents/modelling'
        modelling_home = root  + '/home/Documents/modelling'
        modelling_local = modelling_home

modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)
from plotting_numerical import plot_redgreen_contrast
# from sendmail import *

from tqdm import tqdm
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################
# %matplotlib inline

def fourier2d(final_concentration):
    # print('Normalised PowerSpectrum of original')
    pixels=final_concentration[5]
    pixels[pixels == 0] = 0.0000000001
    pixels/=np.sum(pixels) #line added
    #fourier transform
    f = np.fft.ifftshift(pixels)
    f = np.fft.fft2(f)
    f = np.fft.fftshift(f)
    if plot==True:
        plt.subplot(142)
        plt.imshow(np.real(f))
        plt.colorbar()
        plt.subplot(143)
        plt.imshow(np.real(f))
        plt.colorbar()
        # max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
    PowerSpectrum = np.real(f)**2 +  np.imag(f)**2
    # PowerSpectrum[max_coordinates]=0
    if plot==True:
        plt.subplot(144)
        plt.imshow(PowerSpectrum)
        plt.colorbar()
    if np.sum(PowerSpectrum)!=0:
        PowerSpectrum/=np.sum(PowerSpectrum)
    # print('sum',np.sum(PowerSpectrum))
    # print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))
    # plt.subplot(132)
    # plt.imshow(PowerSpectrum)#,norm=LogNorm())
    # plt.colorbar()
    # binwidth = 0.00001
    PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
    # print(np.sum(PowerSpectrum_vector))
    # plt.subplot(133)
    # hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
    # plt.yscale('log')
    # plt.show()
    return PowerSpectrum_vector

def entropy(histogram):
    x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
    return x

lhs=True
if lhs==True:
    folder = 'fullcircuit_5716gaussian/1M_turingI'
    variant=0
else:
    var=float(sys.argv[1])
    # var=0.23
    folder = 'fullcircuit_5716gaussian/var%s'%var
    variant='5716gaussian'

circuit_n=2
shape='square'
mechanism = 'fullcircuit'


L=5; x_gridpoints =10; J = L*x_gridpoints


T =2000; t_gridpoints = 10; N = T*t_gridpoints
data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )
parID_ps = {}
plot=False
# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
# for parID in tqdm(parID_list[:20], disable=False):

for parID in tqdm([805,686,472,252,688], disable=True):
    # print('parID',parID)
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
    final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
    if not np.isnan(final_concentration0).any():
        final_concentration = np.round(final_concentration0,4)
        if plot==True:
            plt.figure(figsize=[14,2])
            plt.subplot(141)
            plt.imshow(final_concentration[5])
            plt.colorbar()
        hist_PowerSpectrum  = fourier2d(final_concentration)
        ps = entropy(hist_PowerSpectrum)
        parID_ps[parID] = ps
    # print(final_concentration0)

if lhs==True:
    filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
    filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)

pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )
# if root=='/rds/general/user/mo2016':
#     sendmail('general_filename')

/rds/general/user/mo2016

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
~/Documents/modelling/3954/numerical_confocal/code/entropy/plot_all_colony.py in <module>
    110     else:
    111         filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)
--> 112     final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
    113     if not np.isnan(final_concentration0).any():
    114         final_concentration = np.round(final_concentration0,4)

FileNotFoundError: [Errno 2] No such file or directory: '/rds/general/user/mo2016/home/Documents/modelling/3954/numerical_confocal/results/simulation/square/fullcircuit_5716gaussian/1M_turingI/2Dfinal_circuit2_variant0_square_fullcircuitID805_L5_J50_T2000_N20000.pkl'

#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)

if root == '/Users/mo2016':
    modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
    modelling_home = '/Volumes/mo2016/home/Documents/modelling'
    modelling_local = root + '/Documents/modelling'
    import matplotlib as mpl
    mpl.use('tkagg')

if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
        modelling_ephemeral = root + '/ephemeral/Documents/modelling'
        modelling_home = root  + '/home/Documents/modelling'
        modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)


from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm

import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################

%matplotlib inline

# %matplotlib inline
def fourier2d(final_concentration):
    # print('Normalised PowerSpectrum of original')
    pixels=final_concentration[5]
    pixels[pixels == 0] = 0.0000000001
    pixels/=np.sum(pixels) #line added



    #fourier transform
    f = np.fft.ifftshift(pixels)
    f = np.fft.fft2(f)
    f = np.fft.fftshift(f)
    if plot==True:
        plt.subplot(142)
        plt.imshow(np.real(f))
        plt.colorbar()
        plt.subplot(143)
        plt.imshow(np.real(f))
        plt.colorbar()
        # max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
    PowerSpectrum = np.real(f)**2 +  np.imag(f)**2
    # PowerSpectrum[max_coordinates]=0
    if plot==True:
        plt.subplot(144)
        plt.imshow(PowerSpectrum)
        plt.colorbar()

    if np.sum(PowerSpectrum)!=0:
        PowerSpectrum/=np.sum(PowerSpectrum)

    # print('sum',np.sum(PowerSpectrum))

    # print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))

    # plt.subplot(132)
    # plt.imshow(PowerSpectrum)#,norm=LogNorm())
    # plt.colorbar()


    # binwidth = 0.00001
    PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
    # print(np.sum(PowerSpectrum_vector))
    # plt.subplot(133)
    # hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
    # plt.yscale('log')
    # plt.show()



    return PowerSpectrum_vector

def entropy(histogram):
    x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
    return x

lhs=True

if lhs==True:
    folder = 'fullcircuit_5716gaussian/1M_turingI'
    variant=0
else:
    var=float(sys.argv[1])
    # var=0.23
    folder = 'fullcircuit_5716gaussian/var%s'%var
    variant='5716gaussian'

circuit_n=2
shape='square'
mechanism = 'fullcircuit'

L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints

data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )

parID_ps = {}

plot=False
# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
    # print('parID',parID)
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)

    final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
    if not np.isnan(final_concentration0).any():
        final_concentration = np.round(final_concentration0,4)
        if plot==True:
            plt.figure(figsize=[14,2])
            plt.subplot(141)
            plt.imshow(final_concentration[5])
            plt.colorbar()
        hist_PowerSpectrum  = fourier2d(final_concentration)
        ps = entropy(hist_PowerSpectrum)

        parID_ps[parID] = ps
    # print(final_concentration0)


if lhs==True:
    filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
    filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)

pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )


# if root=='/rds/general/user/mo2016':
#     sendmail('general_filename')

 15%|█▌        | 3/20 [00:00<00:00, 23.29it/s]

/rds/general/user/mo2016

100%|██████████| 20/20 [00:00<00:00, 22.90it/s]

#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)

if root == '/Users/mo2016':
    modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
    modelling_home = '/Volumes/mo2016/home/Documents/modelling'
    modelling_local = root + '/Documents/modelling'
    import matplotlib as mpl
    mpl.use('tkagg')

if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
        modelling_ephemeral = root + '/ephemeral/Documents/modelling'
        modelling_home = root  + '/home/Documents/modelling'
        modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)


from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm

import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################

%matplotlib inline
plot=True

# %matplotlib inline
def fourier2d(final_concentration):
    # print('Normalised PowerSpectrum of original')
    pixels=final_concentration[5]
    pixels[pixels == 0] = 0.0000000001
    pixels/=np.sum(pixels) #line added



    #fourier transform
    f = np.fft.ifftshift(pixels)
    f = np.fft.fft2(f)
    f = np.fft.fftshift(f)
    if plot==True:
        plt.subplot(142)
        plt.imshow(np.real(f))
        plt.colorbar()
        plt.subplot(143)
        plt.imshow(np.real(f))
        plt.colorbar()
        # max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
    PowerSpectrum = np.real(f)**2 +  np.imag(f)**2
    # PowerSpectrum[max_coordinates]=0
    if plot==True:
        plt.subplot(144)
        plt.imshow(PowerSpectrum)
        plt.colorbar()

    if np.sum(PowerSpectrum)!=0:
        PowerSpectrum/=np.sum(PowerSpectrum)

    # print('sum',np.sum(PowerSpectrum))

    # print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))

    # plt.subplot(132)
    # plt.imshow(PowerSpectrum)#,norm=LogNorm())
    # plt.colorbar()


    # binwidth = 0.00001
    PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
    # print(np.sum(PowerSpectrum_vector))
    # plt.subplot(133)
    # hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
    # plt.yscale('log')
    # plt.show()



    return PowerSpectrum_vector

def entropy(histogram):
    x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
    return x

lhs=True

if lhs==True:
    folder = 'fullcircuit_5716gaussian/1M_turingI'
    variant=0
else:
    var=float(sys.argv[1])
    # var=0.23
    folder = 'fullcircuit_5716gaussian/var%s'%var
    variant='5716gaussian'

circuit_n=2
shape='square'
mechanism = 'fullcircuit'

L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints

data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )

parID_ps = {}

# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
    # print('parID',parID)
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)

    final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
    if not np.isnan(final_concentration0).any():
        final_concentration = np.round(final_concentration0,4)
        if plot==True:
            plt.figure(figsize=[14,2])
            plt.subplot(141)
            plt.imshow(final_concentration[5])
            plt.colorbar()
        hist_PowerSpectrum  = fourier2d(final_concentration)
        ps = entropy(hist_PowerSpectrum)

        parID_ps[parID] = ps
    # print(final_concentration0)


if lhs==True:
    filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
else:
    filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)

pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )


# if root=='/rds/general/user/mo2016':
#     sendmail('general_filename')

  5%|▌         | 1/20 [00:00<00:02,  7.11it/s]

/rds/general/user/mo2016

100%|██████████| 20/20 [00:02<00:00,  7.06it/s]

#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)

if root == '/Users/mo2016':
    modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
    modelling_home = '/Volumes/mo2016/home/Documents/modelling'
    modelling_local = root + '/Documents/modelling'
    import matplotlib as mpl
    mpl.use('tkagg')

if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
        modelling_ephemeral = root + '/ephemeral/Documents/modelling'
        modelling_home = root  + '/home/Documents/modelling'
        modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)


from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm

import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################

%matplotlib inline
plot=True

# %matplotlib inline
def fourier2d(final_concentration):
    # print('Normalised PowerSpectrum of original')
    pixels=final_concentration[5]
    pixels[pixels == 0] = 0.0000000001
    pixels/=np.sum(pixels) #line added



    #fourier transform
    f = np.fft.ifftshift(pixels)
    f = np.fft.fft2(f)
    f = np.fft.fftshift(f)
    if plot==True:
        plt.subplot(142)
        plt.imshow(np.real(f))
        plt.colorbar()
        plt.subplot(143)
        plt.imshow(np.real(f))
        plt.colorbar()
        max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
    PowerSpectrum = np.real(f)**2 +  np.imag(f)**2
    PowerSpectrum[max_coordinates]=0
    if plot==True:
        plt.subplot(144)
        plt.imshow(PowerSpectrum)
        plt.colorbar()

    if np.sum(PowerSpectrum)!=0:
        PowerSpectrum/=np.sum(PowerSpectrum)

    # print('sum',np.sum(PowerSpectrum))

    # print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))

    # plt.subplot(132)
    # plt.imshow(PowerSpectrum)#,norm=LogNorm())
    # plt.colorbar()


    # binwidth = 0.00001
    PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
    # print(np.sum(PowerSpectrum_vector))
    # plt.subplot(133)
    # hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
    # plt.yscale('log')
    # plt.show()



    return PowerSpectrum_vector

def entropy(histogram):
    x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
    return x

lhs=True

if lhs==True:
    folder = 'fullcircuit_5716gaussian/1M_turingI'
    variant=0
else:
    var=float(sys.argv[1])
    # var=0.23
    folder = 'fullcircuit_5716gaussian/var%s'%var
    variant='5716gaussian'

circuit_n=2
shape='square'
mechanism = 'fullcircuit'

L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints

data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )

parID_ps = {}

# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
    # print('parID',parID)
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)

    final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
    if not np.isnan(final_concentration0).any():
        final_concentration = np.round(final_concentration0,4)
        if plot==True:
            plt.figure(figsize=[14,2])
            plt.subplot(141)
            plt.imshow(final_concentration[5])
            plt.colorbar()
        hist_PowerSpectrum  = fourier2d(final_concentration)
        ps = entropy(hist_PowerSpectrum)

        parID_ps[parID] = ps
    # print(final_concentration0)


# if lhs==True:
#     filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
# else:
#     filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)

pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )


# if root=='/rds/general/user/mo2016':
#     sendmail('general_filename')

  5%|▌         | 1/20 [00:00<00:01,  9.63it/s]

/rds/general/user/mo2016

100%|██████████| 20/20 [00:02<00:00,  8.28it/s]

#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)

if root == '/Users/mo2016':
    modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
    modelling_home = '/Volumes/mo2016/home/Documents/modelling'
    modelling_local = root + '/Documents/modelling'
    import matplotlib as mpl
    mpl.use('tkagg')

if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
        modelling_ephemeral = root + '/ephemeral/Documents/modelling'
        modelling_home = root  + '/home/Documents/modelling'
        modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)


from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm

import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################

%matplotlib inline
plot=True
save=True
if plot==True:
    save==False
# %matplotlib inline
def fourier2d(final_concentration):
    # print('Normalised PowerSpectrum of original')
    pixels=final_concentration[5]
    pixels[pixels == 0] = 0.0000000001
    pixels/=np.sum(pixels) #line added



    #fourier transform
    f = np.fft.ifftshift(pixels)
    f = np.fft.fft2(f)
    f = np.fft.fftshift(f)
    if plot==True:
        plt.subplot(142)
        plt.imshow(np.real(f))
        plt.colorbar()
        plt.subplot(143)
        plt.imshow(np.real(f))
        plt.colorbar()
        max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
    PowerSpectrum = np.real(f)**2 +  np.imag(f)**2
    PowerSpectrum[max_coordinates]=0
    if plot==True:
        plt.subplot(144)
        plt.imshow(PowerSpectrum)
        plt.colorbar()

    if np.sum(PowerSpectrum)!=0:
        PowerSpectrum/=np.sum(PowerSpectrum)

    # print('sum',np.sum(PowerSpectrum))

    # print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))

    # plt.subplot(132)
    # plt.imshow(PowerSpectrum)#,norm=LogNorm())
    # plt.colorbar()


    # binwidth = 0.00001
    PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
    # print(np.sum(PowerSpectrum_vector))
    # plt.subplot(133)
    # hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
    # plt.yscale('log')
    # plt.show()



    return PowerSpectrum_vector

def entropy(histogram):
    x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
    return x

lhs=True

if lhs==True:
    folder = 'fullcircuit_5716gaussian/1M_turingI'
    variant=0
else:
    var=float(sys.argv[1])
    # var=0.23
    folder = 'fullcircuit_5716gaussian/var%s'%var
    variant='5716gaussian'

circuit_n=2
shape='square'
mechanism = 'fullcircuit'

L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints

data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )

parID_ps = {}

# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
    # print('parID',parID)
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)

    final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
    if not np.isnan(final_concentration0).any():
        final_concentration = np.round(final_concentration0,4)
        if plot==True:
            plt.figure(figsize=[24,4])
            plt.subplot(141)
            plt.imshow(final_concentration[5])
            plt.colorbar()
        hist_PowerSpectrum  = fourier2d(final_concentration)
        ps = entropy(hist_PowerSpectrum)

        parID_ps[parID] = ps
    # print(final_concentration0)

if save==True:
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)

pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )


# if root=='/rds/general/user/mo2016':
#     sendmail('general_filename')

  5%|▌         | 1/20 [00:00<00:02,  7.31it/s]

/rds/general/user/mo2016

100%|██████████| 20/20 [00:02<00:00,  7.77it/s]

#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)

if root == '/Users/mo2016':
    modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
    modelling_home = '/Volumes/mo2016/home/Documents/modelling'
    modelling_local = root + '/Documents/modelling'
    import matplotlib as mpl
    mpl.use('tkagg')

if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
        modelling_ephemeral = root + '/ephemeral/Documents/modelling'
        modelling_home = root  + '/home/Documents/modelling'
        modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)


from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm

import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################

%matplotlib inline
plot=True
save=True
if plot==True:
    save==False
# %matplotlib inline
def fourier2d(final_concentration):
    # print('Normalised PowerSpectrum of original')
    pixels=final_concentration[5]
    pixels[pixels == 0] = 0.0000000001
    pixels/=np.sum(pixels) #line added



    #fourier transform
    f = np.fft.ifftshift(pixels)
    f = np.fft.fft2(f)
    f = np.fft.fftshift(f)
    if plot==True:
        plt.subplot(142)
        plt.imshow(np.real(f))
        plt.colorbar()
        plt.subplot(143)
        plt.imshow(np.imag(f))
        plt.colorbar()
        max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
    PowerSpectrum = np.real(f)**2 +  np.imag(f)**2
    PowerSpectrum[max_coordinates]=0
    if plot==True:
        plt.subplot(144)
        plt.imshow(PowerSpectrum)
        plt.colorbar()

    if np.sum(PowerSpectrum)!=0:
        PowerSpectrum/=np.sum(PowerSpectrum)

    # print('sum',np.sum(PowerSpectrum))

    # print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))

    # plt.subplot(132)
    # plt.imshow(PowerSpectrum)#,norm=LogNorm())
    # plt.colorbar()


    # binwidth = 0.00001
    PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
    # print(np.sum(PowerSpectrum_vector))
    # plt.subplot(133)
    # hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
    # plt.yscale('log')
    # plt.show()



    return PowerSpectrum_vector

def entropy(histogram):
    x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
    return x

lhs=True

if lhs==True:
    folder = 'fullcircuit_5716gaussian/1M_turingI'
    variant=0
else:
    var=float(sys.argv[1])
    # var=0.23
    folder = 'fullcircuit_5716gaussian/var%s'%var
    variant='5716gaussian'

circuit_n=2
shape='square'
mechanism = 'fullcircuit'

L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints

data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )

parID_ps = {}

# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
    # print('parID',parID)
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)

    final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
    if not np.isnan(final_concentration0).any():
        final_concentration = np.round(final_concentration0,4)
        if plot==True:
            plt.figure(figsize=[24,4])
            plt.subplot(141)
            plt.imshow(final_concentration[5])
            plt.colorbar()
        hist_PowerSpectrum  = fourier2d(final_concentration)
        ps = entropy(hist_PowerSpectrum)

        parID_ps[parID] = ps
    # print(final_concentration0)

if save==True:
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)

pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )


# if root=='/rds/general/user/mo2016':
#     sendmail('general_filename')

  0%|          | 0/20 [00:00<?, ?it/s]

/rds/general/user/mo2016

100%|██████████| 20/20 [00:03<00:00,  6.41it/s]

#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)

if root == '/Users/mo2016':
    modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
    modelling_home = '/Volumes/mo2016/home/Documents/modelling'
    modelling_local = root + '/Documents/modelling'
    import matplotlib as mpl
    mpl.use('tkagg')

if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
        modelling_ephemeral = root + '/ephemeral/Documents/modelling'
        modelling_home = root  + '/home/Documents/modelling'
        modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)


from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm

import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################

%matplotlib inline
plot=True
save=True
if plot==True:
    save==False
# %matplotlib inline
def fourier2d(final_concentration):
    # print('Normalised PowerSpectrum of original')
    pixels=final_concentration[5]
    pixels[pixels == 0] = 0.0000000001
    pixels/=np.sum(pixels) #line added



    #fourier transform
    f = np.fft.ifftshift(pixels)
    f = np.fft.fft2(f)
    f = np.fft.fftshift(f)
    if plot==True:
        plt.subplot(142)
        plt.imshow(np.real(f))
        plt.colorbar()
        plt.subplot(143)
        plt.imshow(np.imag(f))
        plt.colorbar()
        # max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
    PowerSpectrum = np.real(f)**2 +  np.imag(f)**2
    # PowerSpectrum[max_coordinates]=0
    if plot==True:
        plt.subplot(144)
        plt.imshow(PowerSpectrum)
        plt.colorbar()

    if np.sum(PowerSpectrum)!=0:
        PowerSpectrum/=np.sum(PowerSpectrum)

    # print('sum',np.sum(PowerSpectrum))

    # print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))

    # plt.subplot(132)
    # plt.imshow(PowerSpectrum)#,norm=LogNorm())
    # plt.colorbar()


    # binwidth = 0.00001
    PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
    # print(np.sum(PowerSpectrum_vector))
    # plt.subplot(133)
    # hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
    # plt.yscale('log')
    # plt.show()



    return PowerSpectrum_vector

def entropy(histogram):
    x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
    return x

lhs=True

if lhs==True:
    folder = 'fullcircuit_5716gaussian/1M_turingI'
    variant=0
else:
    var=float(sys.argv[1])
    # var=0.23
    folder = 'fullcircuit_5716gaussian/var%s'%var
    variant='5716gaussian'

circuit_n=2
shape='square'
mechanism = 'fullcircuit'

L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints

data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )

parID_ps = {}

# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
    # print('parID',parID)
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)

    final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
    if not np.isnan(final_concentration0).any():
        final_concentration = np.round(final_concentration0,4)
        if plot==True:
            plt.figure(figsize=[24,4])
            plt.subplot(141)
            plt.imshow(final_concentration[5])
            plt.colorbar()
        hist_PowerSpectrum  = fourier2d(final_concentration)
        ps = entropy(hist_PowerSpectrum)

        parID_ps[parID] = ps
    # print(final_concentration0)

if save==True:
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)

pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )


# if root=='/rds/general/user/mo2016':
#     sendmail('general_filename')

  5%|▌         | 1/20 [00:00<00:02,  7.40it/s]

/rds/general/user/mo2016

100%|██████████| 20/20 [00:02<00:00,  8.44it/s]

#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)

if root == '/Users/mo2016':
    modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
    modelling_home = '/Volumes/mo2016/home/Documents/modelling'
    modelling_local = root + '/Documents/modelling'
    import matplotlib as mpl
    mpl.use('tkagg')

if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
        modelling_ephemeral = root + '/ephemeral/Documents/modelling'
        modelling_home = root  + '/home/Documents/modelling'
        modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)


from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm

import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################

%matplotlib inline
plot=True
save=True
if plot==True:
    save==False
# %matplotlib inline
def fourier2d(final_concentration):
    # print('Normalised PowerSpectrum of original')
    pixels=final_concentration[5]
    pixels[pixels == 0] = 0.0000000001
    pixels/=np.sum(pixels) #line added



    #fourier transform
    f = np.fft.ifftshift(pixels)
    f = np.fft.fft2(f)
    f = np.fft.fftshift(f)
    if plot==True:
        plt.subplot(142)
        plt.imshow(np.real(f))
        plt.colorbar()
        plt.subplot(143)
        plt.imshow(np.imag(f))
        plt.colorbar()
        # max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
    PowerSpectrum = np.real(f)**2 +  np.imag(f)**2
    # PowerSpectrum[max_coordinates]=0
    # if plot==True:
    #     plt.subplot(144)
    #     plt.imshow(PowerSpectrum)
    #     plt.colorbar()

    if np.sum(PowerSpectrum)!=0:
        PowerSpectrum/=np.sum(PowerSpectrum)
    if plot==True:
        plt.subplot(144)
        plt.imshow(PowerSpectrum)
        plt.colorbar()
    # print('sum',np.sum(PowerSpectrum))

    # print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))

    # plt.subplot(132)
    # plt.imshow(PowerSpectrum)#,norm=LogNorm())
    # plt.colorbar()


    # binwidth = 0.00001
    PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
    # print(np.sum(PowerSpectrum_vector))
    # plt.subplot(133)
    # hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
    # plt.yscale('log')
    # plt.show()



    return PowerSpectrum_vector

def entropy(histogram):
    x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
    return x

lhs=True

if lhs==True:
    folder = 'fullcircuit_5716gaussian/1M_turingI'
    variant=0
else:
    var=float(sys.argv[1])
    # var=0.23
    folder = 'fullcircuit_5716gaussian/var%s'%var
    variant='5716gaussian'

circuit_n=2
shape='square'
mechanism = 'fullcircuit'

L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints

data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )

parID_ps = {}

# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
    # print('parID',parID)
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)

    final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
    if not np.isnan(final_concentration0).any():
        final_concentration = np.round(final_concentration0,4)
        if plot==True:
            plt.figure(figsize=[24,4])
            plt.subplot(141)
            plt.imshow(final_concentration[5])
            plt.colorbar()
        hist_PowerSpectrum  = fourier2d(final_concentration)
        ps = entropy(hist_PowerSpectrum)

        parID_ps[parID] = ps
    # print(final_concentration0)

if save==True:
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)

pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )


# if root=='/rds/general/user/mo2016':
#     sendmail('general_filename')

  5%|▌         | 1/20 [00:00<00:01,  9.92it/s]

/rds/general/user/mo2016

100%|██████████| 20/20 [00:02<00:00,  7.57it/s]

#############################
#IMPORTS#
#############################
import sys
import os
pwd = os.getcwd()
root = pwd.rpartition("mo2016")[0] + pwd.rpartition("mo2016")[1] #/Volumes/mo2016/ or '/Users/mo2016/' or '/rds/general/mo2016/'
print(root)

if root == '/Users/mo2016':
    modelling_ephemeral = '/Volumes/mo2016/ephemeral/Documents/modelling'
    modelling_home = '/Volumes/mo2016/home/Documents/modelling'
    modelling_local = root + '/Documents/modelling'
    import matplotlib as mpl
    mpl.use('tkagg')

if root == '/Volumes/mo2016' or root=='/rds/general/user/mo2016': #'/rds/general' or root=='/Volumes':
        modelling_ephemeral = root + '/ephemeral/Documents/modelling'
        modelling_home = root  + '/home/Documents/modelling'
        modelling_local = modelling_home
modulepath = modelling_local + '/3954/modules/new_CN'
sys.path.append(modulepath)


from plotting_numerical import plot_redgreen_contrast
# from sendmail import *
from tqdm import tqdm

import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

from scipy.fft import fft, ifft
import matplotlib as mpl
from PIL import Image, ImageDraw
import numpy as np
from numpy import asarray
import pickle
from tqdm import tqdm
#############################

%matplotlib inline
plot=True
save=True
if plot==True:
    save==False
# %matplotlib inline
def fourier2d(final_concentration):
    # print('Normalised PowerSpectrum of original')
    pixels=final_concentration[5]
    pixels[pixels == 0] = 0.0000000001
    pixels/=np.sum(pixels) #line added



    #fourier transform
    f = np.fft.ifftshift(pixels)
    f = np.fft.fft2(f)
    f = np.fft.fftshift(f)
    if plot==True:
        plt.subplot(142)
        plt.imshow(np.real(f))
        plt.colorbar()
        plt.subplot(143)
        plt.imshow(np.imag(f))
        plt.colorbar()
        max_coordinates=list(zip(np.where(f == np.amax(f))[0], np.where(f == np.amax(f))[1]))[0]
    PowerSpectrum = np.real(f)**2 +  np.imag(f)**2
    PowerSpectrum[max_coordinates]=0
    # if plot==True:
    #     plt.subplot(144)
    #     plt.imshow(PowerSpectrum)
    #     plt.colorbar()

    if np.sum(PowerSpectrum)!=0:
        PowerSpectrum/=np.sum(PowerSpectrum)
    if plot==True:
        plt.subplot(144)
        plt.imshow(PowerSpectrum)
        plt.colorbar()
    # print('sum',np.sum(PowerSpectrum))

    # print('max', np.amax(PowerSpectrum), 'min', np.amin(PowerSpectrum))

    # plt.subplot(132)
    # plt.imshow(PowerSpectrum)#,norm=LogNorm())
    # plt.colorbar()


    # binwidth = 0.00001
    PowerSpectrum_vector = np.reshape(PowerSpectrum,-1)
    # print(np.sum(PowerSpectrum_vector))
    # plt.subplot(133)
    # hist_PowerSpectrum= plt.hist(PowerSpectrum_vector, bins=np.linspace(np.amin(PowerSpectrum_vector), np.amax(PowerSpectrum_vector) + binwidth, 300))
    # plt.yscale('log')
    # plt.show()



    return PowerSpectrum_vector

def entropy(histogram):
    x= np.sum([-p*np.log2(p) for p in histogram if p!=0])
    return x

lhs=True

if lhs==True:
    folder = 'fullcircuit_5716gaussian/1M_turingI'
    variant=0
else:
    var=float(sys.argv[1])
    # var=0.23
    folder = 'fullcircuit_5716gaussian/var%s'%var
    variant='5716gaussian'

circuit_n=2
shape='square'
mechanism = 'fullcircuit'

L=5; x_gridpoints =10; J = L*x_gridpoints
T =2000; t_gridpoints = 10; N = T*t_gridpoints

data_path = modelling_home + '/3954/numerical_confocal/results/simulation/square/%s'%(folder)
parID_list = pickle.load( open(data_path + '/parID_list_L%r_J%r_T%r_N%r.pkl'%(L,J,T,N), "rb" ) )

parID_ps = {}

# parID_list=[497]
# for parID in tqdm(parID_list, disable=False):
for parID in tqdm(parID_list[:20], disable=False):
# for parID in tqdm([805,686,472,252,688], disable=True):
    # print('parID',parID)
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,int(parID),L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%r_%s_%sID%r_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,int(parID),L,J,T,N)

    final_concentration0 = pickle.load( open(data_path + '/2Dfinal_%s.pkl'%filename, 'rb' ) )
    if not np.isnan(final_concentration0).any():
        final_concentration = np.round(final_concentration0,4)
        if plot==True:
            plt.figure(figsize=[24,4])
            plt.subplot(141)
            plt.imshow(final_concentration[5])
            plt.colorbar()
        hist_PowerSpectrum  = fourier2d(final_concentration)
        ps = entropy(hist_PowerSpectrum)

        parID_ps[parID] = ps
    # print(final_concentration0)

if save==True:
    if lhs==True:
        filename = 'circuit%r_variant%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant, shape,mechanism,L,J,T,N)
    else:
        filename = 'circuit%r_variant%svar%s_%s_%s_L%r_J%r_T%r_N%r'%(circuit_n,variant,var, shape,mechanism,L,J,T,N)

pickle.dump( parID_ps, open( modelling_home + "/3954/numerical_confocal/results/entropy/EntropyDicts/ps_zerokept_dict_%s.pkl"%filename, "wb" ) )


# if root=='/rds/general/user/mo2016':
#     sendmail('general_filename')

  5%|▌         | 1/20 [00:00<00:01,  9.63it/s]

/rds/general/user/mo2016

100%|██████████| 20/20 [00:02<00:00,  7.95it/s]